Deep brain stimulation of the thalamus restores signatures of consciousness in a nonhuman primate model

Loss of consciousness is associated with the disruption of long-range thalamocortical and corticocortical brain communication. We tested the hypothesis that deep brain stimulation (DBS) of central thalamus might restore both arousal and awareness following consciousness loss. We applied anesthesia to suppress consciousness in nonhuman primates. During anesthesia, central thalamic stimulation induced arousal in an on-off manner and increased functional magnetic resonance imaging activity in prefrontal, parietal, and cingulate cortices. Moreover, DBS restored a broad dynamic repertoire of spontaneous resting-state activity, previously described as a signature of consciousness. None of these effects were obtained during the stimulation of a control site in the ventrolateral thalamus. Last, DBS restored a broad hierarchical response to auditory violations that was disrupted under anesthesia. Thus, DBS restored the two dimensions of consciousness, arousal and conscious access, following consciousness loss, paving the way to its therapeutical translation in patients with disorders of consciousness.

protected with a plastic MR compatible chamber that was home-made by 3D-printing and was fixed to the skull with screws and dental acrylic.

Verification of the DBS settings and underlying behavioral responses
To ensure efficiency across experiments, impedances between leads and through electrode to an external reference were first measured outside the MRI environment with the DBS programmer device provided by the manufacturer (8840 N'Vision, Activa Clinician DBS Programmer, Medtronic, USA).
Ultimately, we used an oscilloscope (Wave Runner 44XI, LeCroy, USA) to check the electrical current delivered to each lead at the beginning and at the end of each experimental session. The stimulation also generated an artifact on the EEG signal, which provided a final benchmark during the fMRI acquisitions.
Behavioral responses (see behavioral assessment section) to electrical thalamic stimulation were assessed in each animal outside the scanner at least 20 days after the DBS implantation. We empirically explored different voltage amplitudes and pulse widths while keeping a monopolar stimulation at a frequency of 130Hz, applied successively to each of the four DBS contacts. The DBS lead targeted the centro-median thalamus. On the target contact (centered in CM), we determined the voltage level for high central thalamic (CT) DBS as the voltage just above the threshold at which a significant behavioral response was observed. Low CT-DBS corresponded to a lower current delivery below the voltage level that led to an arousal pattern. For comparison and reproducibility purposes, we kept the exact same DBS settings for the control stimulation site (ventro-lateral thalamus (VL) DBS).

Block-design fMRI analysis
We generated plots by extracting the activations responses to high CT-DBS with the hemodynamic function in frontal (area 6V; 9/46; 8A), parietal (ventral intraparietal VIP, parietal area PFG), cingulate (anterior ACC: posterior PCC) and temporal cortex (temporo-parieto-occipital area TPO). Activity profiles were plotted as percentage of signal changes across time.

Resting-state fMRI analysis
For the static resting state, we calculated for each experimental condition and sessions the average of positive and negative Z-values and performed a Student t-test with the null hypothesis of zero correlation to test for statistical significance of connectivity between the different experimental conditions.
We computed the static functional correlations by estimating for each experimental condition (noted e) for the awake state, anesthesia, low CT-DBS, high CT-DBS, low VL-DBS and high VL-DBS and acquired run (noted r) the covariance matrix Ce,r. This value was obtained by extracting and averaging across all runs r the time series of all voxels included in each selected anatomical ROI. We referred as static functional correlation or stationary functional connectivity the entry matrix Ce,r (i,j) where each cell represented the mean strength of the functional correlation between the i-j pair. For each covariance matrix, a Fisher transformation was applied to calculate the Z-score. The Z-score matrices were averaged across runs to obtain one matrix per experimental condition. To assess statistical significance, Student t-test at the threshold p value<0.0001 and a false discovery rate correction were applied on the correlations of all pairs of brain regions.

Dynamic resting state fMRI analysis
Covariance values between all ROIs were included ([82×(82−1)]/2=3,321 features per matrix). Zc,s,w matrices were subsampled along the time dimension (w) before clustering. The resulting centroids or median clusters (BSn with n=1-7; each BSn is sized 82×82) were then used to initialize a clustering of all data, obtaining a matrix of brain states Bc,s,w, which, for a given arousal condition c and session s, is a vector of length 464, valued 1-7, because each matrix in Zs,p,w is assigned a BSn.
The similarity score was computed from the correlation coefficient between the vectorized structural matrix and each vectorized brain state from the clustering analysis. All brain states were ordered in ascending order of similarity to the structure using the similarity score. To quantify the relation between the probability of occurrence of a BS and the similarity score, for each arousal condition, a regression analysis was done, to quantify the beta value (β), the R 2 and a P value. The differences in BS composition across arousal states was evaluated through a fixed-effects ANOVA, with mean rank similarity, that is, the result of averaging each BS time series, valued from 1 to 7, as a dependent variable and the vigilance condition as the in-dependent variable. A fixed-effects ANOVA was run to quantify the effect of sedation on the probability of brain state 7. For this, we followed the same procedure, but the mean rank similarity was calculated considering only BS 7 (window w valued at 1; any other state window w valued at 0).
To explore specifically the fluctuations of intervoxel correlation within nodes of the "macaque global neuronal worskspace" and sensori-motor areas (anterior cingulate cortex, ACC; dorsolateral prefrontal cortex, PFCdl; frontal eye fields, FEF; dorsolateral premotor cortex, PMCdl; primary somatosensory cortex, S1; primary motor cortex, M1; intraparietal cortex, Pcip; primary auditory cortex, A1; inferior temporal cortex, TCi; visual area 1, V1 and posterior cingulate cortex, PCC), we extracted the values from the whole brain matrices and applied a one-way analysis of variance (ANOVA). FC across the brain states were highlighted by displaying Z-score in inter-region matrices.

Event-related task fMRI analysis
We generated plots by obtaining the β-weight of SPM regressions of individual macaque data with the hemodynamic functions of the appropriate stimulus categories and then plotted the mean and SE of these β-weights. These values estimate, in percentages of the whole-brain fMRI signal, the size of the fMRI activation relative to the implicit rest baseline that divides trials.
The first level analyses consisted in the convolution of the stimulus categories with the MION canonical hemodynamic response function (HRF) and its time derivative. We also added motion regressors and heart rate as variables of non-interest to the event-related regressors. Activation time series of all the fMRI voxels were computed for each fMRI run and signal change expressed in Tscore maps for the different stimulus categories relative to rest periods. Global standard trials that immediately followed a global deviant trial were excluded.

-Extended Results
Thalamic DBS effects on resting-state networks in anesthetized macaques

Static Functional Correlations (Figure 4, S6)
To test for statistical significance of connectivity between brain regions in different experimental conditions, Student t-tests were performed with the null hypothesis of zero correlation. We calculated for each experimental condition c and sessions the average of positive and negative Z values of Zc,s.

Dynamic Functional Correlations (Figure 5, S7-S8, Table S4)
We applied k-means to the whole acquired dataset (including all experimental conditions) to cluster brain states ( Figure S7). We also applied k-means to two data subsets, subset CT and subset VL, to specifically characterize the effects of CT-DBS and VL-DBS respectively. Subset CT included data from awake, anesthesia and anesthesia + high CT-DBS conditions ( Figure 5A-C). Subset VL included data from awake, anesthesia and anesthesia + high VL-DBS conditions ( Figure 5D-E).
-Clustering subset CT (data from awake, anesthesia and anesthesia + high CT-DBS conditions)  In the awake state, all 7 brain states were represented with a similar probability of occurrence (β=0.45; R 2 =0.22; p=0.28). During anesthesia, state 7 (with the highest function-structure similarity) was dominant and state 1 (with the lowest function-structure similarity) never occurred (β=1.91; R 2 =0.67; p=0.02). Under high CT-DBS, consistent with partial recovery of consciousness, the probability of occurrence of state 7 decreased in favor of all the other brain states, especially state 2 and 3 (β=0.64; R 2 =0.28; p=0.22). We computed the slope of the linear relation between structural and functional correlations for each recording session and compared the slope distributions in the awake, anesthesia and DBS conditions. Awake and high CT-DBS slopes were significantly lower than the anesthesia slopes, indicating that a greater diversity of states were explored in the wake state (awake versus anesthesia: t-test, p=6e-5 and BF10=338, high CT-DBS versus anesthesia: t-test, p=0.001 and BF10=23). Importantly no differences were observed between awake and high CT-DBS slopes (t-test, p=0.42, BF01=3.28) ( Figure 5A-B, Table S6). In the awake state, the mean rank of brain states was 4 (4.38±1.28), for anesthesia, the mean rank was 6 (5.70±1.40) and during high CT-DBS, the mean rank was 5 (4.55±1.17). This brain state distribution was significantly different (ANOVA; F(2;120)=12.52; p=1.15e-7). Also, the frequency of brain state 7 was moderate in the awake experiments (probability=0.24), high during anesthesia (probability=0.58) and low again during high CT-DBS (probability=0.26; p<0.0001) ( Figure 5B). The probability of brain state 7 was higher in anesthesia compared to the awake state (t-test, p=1e-6, BF10=9063) and to high CT-DBS (t-test, p=1e-5, BF10=1017) which did not differ significantly (t-test, p=0.74, BF01=4.18). The mean similarity with the anatomical connectivity was also significantly different with 0.24 (±0.06) for the awake state, 0.31 (±0.07) for anesthesia and 0.25 (±0.06) for high CT-DBS (ANOVA; F(2,120)=14.75; p=1.87e-6).
Anatomically, the functional brain states 1, 2 and 3, that were most characteristic of the awake, presented strong correlations within the "macaque GNW" prefrontal (dorsolateral prefrontal cortex, PFCdl; dorsolateral premotor cortex, PMCdl), parietal (intraparietal cortex, PCip) and cingulate nodes (anterior cingulate cortex, ACC; posterior cingulate cortex, PCCr), whereas state 7 displayed low or null Z-score values across the same entire cortical network ( Figure 5C). During high CT-DBS, the average duration of brain state 7 decreased compared to anesthesia (high CT-DBS versus anesthesia, p=2.48e-3; bootstrap analysis) and was similar to the awake state ( Figure S8).
With high VL-DBS, the average duration of brain state 7 increased compared to the awake state (high VL-DBS v/s awake, p<0.0001, bootstrap analysis) and was similar to the anesthesia state ( Figure S8).
-Clustering the whole dataset (data from awake, anesthesia, low CT-DBS, high CT-DBS, low

C.
A. B.

D.
A.

B.
Right ACC Right PFC 6V F.  B.

(CT) DBS (left) or high ventral-lateral thalamic (VL) DBS (right); normalized probability distribution and twodimensional normalized histograms
Average life time of brain states for the awake, anesthesia and high central thalamic (CT) DBS condition (A) or high ventrallateral thalamic (VL) DBS (B) for 7 brain states obtained by k-means clustering. Error bars stand for 1 SEM. Normalized probability distribution of all Z values for the functional brain state 1 (the least similar to the structural brain connectivity) and functional brain state 7 (the most similar to the structural brain connectivity) for awake, anesthesia and high CT-DBS resting state pooled together. Similar results were obtained regardless the inputs conditions for the clustering (C). Two-dimensional normalized histograms for functional brain state 1 and functional brain state 7 for the clustering of awake, anesthesia and high CT-DBS condition. (D) Z values as a function of distance between pairs of regions of interest for brain state 1 (upper right) and brain state 7 (lower right) for the clustering of awake, anesthesia and high CT-DBS condition.  Local deviants occur at the trial level (1 st order) whereas Global deviants occur at the series level (2 nd order).
A.     A.      Oxygen saturation (SpO 2 ); systolic, diastolic and mean blood pressure (respectively SBP, DBP, MBP); respiration rate (RR); end-tidal CO 2 (EtCO 2 ) and temperature (T) for each animal (monkey N and T) under general anesthesia, general anesthesia plus low central thalamic (CT) DBS, general anesthesia plus high CT-DBS, general anesthesia plus low ventral-lateral thalamic (VL) DBS and general anesthesia plus high VL-DBS.  Value of the Bayes Factor BF10 and BF01 to interpret statistical evidence in favor of the H1 or H0 hypothesis. A BF greater than 3 significatively support the evidence of the tested hypothesis.